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Abstract 

This paper investigates the state estimation of a high-hdelity spatially resolved thermal- 
electrochemical lithium-ion battery model commonly referred to as the pseudo two-dimensional 
model. The partial-differential algebraic equations (PDAEs) constituting the model are spa¬ 
tially discretised using Chebyshev orthogonal collocation enabling fast and accurate simu¬ 
lations up to high C-rates. This implementation of the pseudo-2D model is then used in 
combination with an extended Kalman hlter algorithm for differential-algebraic equations to 
estimate the states of the model. The state estimation algorithm is able to rapidly recover 
the model states from current, voltage and temperature measurements. Results show that 
the error on the state estimate falls below 1 % in less than 200 s despite a 30 % error on bat¬ 
tery initial state-of-charge and additive measurement noise with 10 mV and 0.5 K standard 
deviations. 

Keywords: Lithium-ion battery, pseudo-two dimensional model, state estimation, 
extended Kalman hlter, Chebyshev orthogonal collocation 


1. Introduction 

Lithium-ion batteries are widely used in electric vehicles and hybrid electric vehicles due 
their high energy and power density compared to other battery chemistries, and are increas¬ 
ingly of interest in grid and oh-grid applications. However, scaling up the size of battery 
packs for automotive and other applications raises new safety and reliability challenges that 
require development of novel sophisticated battery management systems (BMSs). A BMS 
consists of hardware and embedded algorithms that ensure the safe and reliable operation 
of a pack by monitoring cells and estimating their states, such as state-of-charge (SOC) 
and state-of-health (SOH) [1]. In order to infer unmeasurable states from the available 
measurements of voltage, current and temperature, a model must be solved in the BMS. 
In automotive applications, the model should accurately describe behaviour under the wide 
range of operating conditions encountered, including high current, extreme temperatures 
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and highly dynamic loads. In addition, diagnosis and prognosis of degradation in terms of 
capacity and power fade is an acute challenge. 

Current BMSs typically employ low-order empirical models, such as equivalent-circuit 
models (ECMs), which are parametrised using time- or frequency-domain experimental data 
[2, 3, 4] for battery state estimation and control. These models have relatively low computa¬ 
tional demands but are only valid within the narrow operating conditions in which they have 
been parametrised. Because the parameters of such models have little physical signihcance, 
broadening their validity range requires a large amount of experimental data under a wide 
range of operating conditions, and predicting degradation is challenging or impossible. 

Alternatively, physics-based models describing the thermodynamics, reaction kinetics and 
transport within the cell are valid over a wide range of operating conditions and could be 
coupled to degradation models directly. Physics-based models have been widely used for 
battery design [5, 6, 7, 8] but are usually too computationally intensive for the limited re¬ 
sources of an embedded BMS. The so-called pseudo two-dimensional (P2D) model developed 
by the Newman group [5] is probably the most widely used lithium-ion battery model of this 
type. It is composed of a one-dimensional macro-scale model describing the evolution of 
lithium concentration and electric potential in the electrolyte across the anode, separator 
and cathode and micro-scale models for the electrodes. The pseudo-second dimension arises 
from these coupled one-dimensional micro-scale models describing the solid-phase diffusion 
of lithium in the porous active material of the electrodes. These micro-scale models solve the 
diffusion of lithium occurring in a spherical particle at each local position of the macro-scale 
porous electrode model. By modelling diffusion and kinetics limitations, the P2D model is 
able to accurately describe lithium-ion battery dynamics over a wide operating range [9] and 
is therefore an excellent starting point for the next generation of BMSs. 

However, the computation required by the P2D model is intense compared to ECMs 
for embedded applications. Several attempts at performing state estimation on simplihed 
models derived from the P2D model have been reported in the literature. A common sim- 
plihcation known as the single particle model (SPM), assumes that each electrode can be 
represented by a unique solid-phase particle and neglects concentration gradients in the 
electrolyte. State estimation using the SPM and similar approximations has already been 
reported in the literature, and includes the use of an extended Kalman filter (EKE) algo¬ 
rithm [10], or a backstopping PDE state estimator [11]. In [12], the EKE was applied to an 
averaged electrochemical model similar to the SPM to estimate SOC and critical concen¬ 
tration at the surface of the electrodes. However, these approaches are inherently limited 
due to the low current validity range of the SPM. Other approaches include state estimation 
on reduced-order models derived from the P2D model. In [13, 14], Kalman hltering is per¬ 
formed on a reduced-order state variable model computed by residue grouping [15, 16] from 
transcendental transfer functions approximating each equation of the P2D model assuming 
quasi-linear behaviour. In [17], the EKE is applied to a state space reduced-order model 
computed from the P2D model using a discrete-time realization algorithm [18, 19, 20]. How¬ 
ever, the parameters of such reduced-order models may be difficult to interpret or have no 
direct physical meaning, which makes accounting for degradation effects difficult. 

Recent works have shown that using spectral numerical methods instead of the commonly 
used hnite-difference method to discretise the P2D model results in a highly reduced model 
order whilst maintaining accuracy and physical signihcance of parameters. Dao et ah used 
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the Galerkin spectral method on sinusoidal basis functions to discretise the electrolyte dif¬ 
fusion equation [21], while Cai and White applied orthogonal collocation on hnite elements 
to all the equations of the P2D model [22], Orthogonal collocation enforced at zeros of 
Jacobi polynomials was also applied to the full P2D model in [23] and solved using Maple 
and DASSL solvers using cosine basis functions and more recently Chebyshev polynomial 
basis functions for improved convergence at high currents [24]. In previous work, we applied 
Chebyshev orthogonal collocation to the isothermal P2D model and showed that compu¬ 
tation time could be reduced by a factor of 10 to 100 compared to hnite-difference for the 
same result accuracy [25]. We have also successfully applied this approach for simulation of 
supercapacitors [26]. 

In this work, we applied the EKF algorithm to the thermal-electrochemical P2D model 
solved using Chebyshev orthogonal collocation for battery state estimation. State estimation 
of the full P2D model solved using the approach discussed in [23] has recently been reported 
using the optimisation-based moving-horizon estimation technique [27] and a tethered par¬ 
ticle hlter algorithm [28]. However, our approach using the simpler EKF algorithm is much 
less computationally intensive while showing good performance. To our knowledge, this is 
the first attempt at estimating the states of the full P2D model using the EKF algorithm. 
Solving a high-fidelity model such as the P2D model coupled to degradation models online 
a BMS can provide valuable information on the internal states of the battery, enabling new 
safety limits [29] and advanced health-conscious control algorithms [30] to be used. 

2. Thermal-electrochemical model 

The battery model considered here consists of the P2D electrochemical model coupled 
to a bulk thermal model. The electrochemical model describes the transport of lithium, 
reaction kinetics and thermodynamics at the electrode level while the bulk thermal model 
describes the evolution of cell temperature. The electrochemical and thermal models are 
coupled together through the potential- and concentration-dependent heat generation rate 
and the temperature-dependent physical and chemical properties of the P2D model. 

2.1. Electrochemical model 

A lithium-ion cell consists of two porous electrodes composed of an active material that 
can store lithium intercalated in the solid material, and a separator that allows the passage 
of ions but not electrons. The electrodes and the separator are soaked in an electrolyte that 
allows the transport of ions. During discharge, lithium stored in the anode is de-inserted from 
the active material and released as ions in the electrolyte. Driven by diffusion (concentration 
gradient) and migration (potential gradient), lithium ions travel through the separator to the 
cathode where they are inserted in the lattice of the cathode active material. Simultaneously, 
electrons travel from the anode to the cathode through the external circuit, powering a load, 
to ensure electro-neutrality. This process is reversed during battery charging. 

From a mathematical modelling perspective, the cell is divided into three domains: anode, 
separator and cathode, denoted Da, D* and Dc respectively (Fig. 1). In each of these domains, 
two phases are considered, the solid phase and the electrolyte phase, and are treated as 
superimposed continua using porous electrode theory [31], sometimes called homogenization, 
therefore neglecting the exact micro-structure of the electrodes. In order to account for the 
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Figure 1: Schematic of the cell computational domains. The cell is divided in three domains, anode fla, 
separator 17^ and cathode 17^, where two phases are present: the electrolyte phase and the solid phase. The 
porous nature of each electrode is considered by assuming spherical particles of solid-phase material fla,p and 
f7c,p at each local position in the anode and cathode domains respectively. The physical coordinate of the 
domain are the j;-coordinates across the cell and the radial r-coordinates in the particle. The computational 
domain is rescaled to [—1,1] and consists of different sets of Chebyshev collocation nodes Xa, Xs, Xc and Xp. 


tortuosity of the porous material, effective electrolyte diffusivity and ionic 

conductivity = K,el are considered with b the so-called Bruggeman coefficient [5, 7, 32]. 

The P2D model consists of a set of partial differential equations (PDEs) and algebraic 
constraints governing the evolution of lithium concentration and electric potential within the 
cell. The dependent variables are solid-phase concentration Cs{r,x,t), electrolyte concentra¬ 
tion Ce{x,t), electric potential at the surface of the solid-phase particles (()s{,x,t), electric 
potential in the electrolyte t) and volumetric reaction current j^^{x, t), which expresses 
the amount of lithium exchanged between the solid-phase and the electrolyte per unit volume 
of electrode. The independent variables are time t, the x-coordinate across the cell thick¬ 
ness and the spherical r-coordinate in the solid-phase particles. The transport of lithium 
in each spherical particle is described by the spherical diffusion equation (1) with Neumann 
boundary conditions (2): 
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The change of variable Cg = rcg is introduced to simplify this sub-model implementation (A). 
The reaction current is given as a function of the electrode local overpotential rj by the 
Butler-Volmer kinetics equation: 




Ugio 


f OlaT \ 

exp I 1 - exp 



(3) 


where = Seg/Rg is the specihc interfacial area of the electrode. The exchange current den¬ 
sity io depends on the particle surface concentration and the electrolyte concentration 
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Ce according to: 


io = kJ^ (cr^)“^ (Ce)““ (4) 

The overpotential in (3) is given hj r] = (ps — 4>e — U where the experimentally-htted 
open-circuit potential functions 17(c®“'’4) are taken from [33]. 

The evolution of lithium concentration in the electrolyte is governed by the diffusion 
equation (5) subject to homogeneous Neumann boundary conditions. 
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The electrolyte potential (pe is governed by Ohm’s law: 
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where the diffusional conductivity is given by: 
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(7) 


ai’ = ^ (1 - 2«"+) an (8) 

It has recently been reported [34] that mistakes are sometimes made in the literature re¬ 
garding (7) and (8) and we therefore are careful to use the correct expressions here. The 
electrolyte potential at the cathode current collector is chosen as the reference potential and 
set to zero to ensure that the system of equations is fully constrained. The solid phase 
potential at the surface of the particles is governed by Ohm’s law: 

-f = 0 (9) 

The local fractions of current density carried by the ions in electrolyte if. and electrons in 
the solid-phase is are related to the total current density passing through the cell iapp by 
the Kirchoff’s law is + ie = iapp- By virtue of conservation of charge, the local reaction rate 
is equal to the divergence of the electrolyte current density. The input of the model 
is the current I from which the applied current density iapp = I/Ag is calculated knowing 
the electrode surface area Ag. The terminal voltage of the cell V is equal to the difference 
between the solid-phase potential at the cathode current collector and that at the anode 
current-collector, minus the ohmic drop due to the contact resistance Rc = 20 kl.crri^ [35] at 
the current collector/electrode interfaces. 

The cell considered for this study consists of a lithium cobalt oxide (LiCo02) cathode and 
a mesocarbon microbead (MCMB) anode with IM LiPFg in propylene carbonate, ethylene 
carbonate and dimethyl carbonate (PC:EC:DMC) electrolyte. The parameters were found in 
the literature and are summarised in Table 1. It has been shown that electrolyte properties 
are highly dependent on lithium concentration and cell temperature [36]. The empirical 
expressions for diffusivity Df. and ionic conductivity k as a function of concentration and 
temperature reported in [36] were used in the present work and the transference number 
Iq = 0.435 is assumed constant [37]. The focus of the present work is on the efficient 
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solution and state estimation of the P2D model. We acknowledge that the estimation of 
model parameters from experimental data is crucial for the practical implementation of a 
state estimation algorithm. This is a challenging task due to the large number of parameters 
compared to the limited number non-invasive measurements available and this will be the 
focus of future work. 


Table 1: Set of parameters of the electrochemical P2D model used for the simulations 


Parameter 

Units 

Anode 

LixCg 

Separator 

LiPFe 

Cathode 

LiyCo02 

Ref. 

5* 

pm 

73.5 

25.0 

70.0 

[22] 

Ri 

pm 

12.5 

— 

8.5 

[22] 

e* 

— 

0.4382 

0.45 

0.3 

[37] 


— 

0.0566 

— 

0.15 

[37] 

Oii 

— 

0.5 

— 

0.5 

[8] 


m?'^ .mol~^'^ .s~^ 

1.764 X 10-1^ 

— 

6.667 X 10“^^ 

[37] 

Ds,i 

2 —1 
m .s 

5.5 X 10“^^ 

— 

1.0 X 10-^1 

[8] 

CTi 

S.m~^ 

100 

— 

10 

[8] 

h 

— 

4.1 

2.3 

1.5 

[37] 

^max 

^5,2 

mol.m~^ 

30555 

— 

51555 

[33] 


— 

0.756 

— 

0.465 

[22] 


2.2. Lumped thermal model 

The P2D electrochemical model is coupled to a lumped thermal model described by the 
following energy balance equation: 


pCp 


dT 

dt 


Qgen T Q 


conv 


( 10 ) 


The total heat generation rate per unit volume g^en is assumed uniform and attributed to 
four main contributions according to g^en = Qrxn + ^rev + Qohm + 9c, where g^^n is the reaction 
heat generation rate, q^ev is the reversible heat generation rate due to entropy changes in the 
active material of electrodes during the intercalation/de-intercalation of lithium and qohm is 
the electronic and ionic ohmic heat generation rate due to the motion of lithium. The average 
heat generated by each of these per unit volume [37, 38] is given by equations (11)-(12)-(13): 


Qrxn 


Qrev 


Qohm 


ij 3^\<t>s-4>e-U°^ndx 


a 


eff 


dx 




dx 

<91nCe \ / ^ 
dx ) \ dx 


dx 


( 11 ) 

( 12 ) 


(13) 


6 



The ohmic heat generated per unit volume due to the contact resistance between the 
electrodes and current collectors is given by: 


Qc = 



(14) 


The rate of convective heat removal per unit volume from the cell to the coolant air qconv in 
(10) is given by: 

_ hAc {T -T^ ) 

Hconv — -rr 

In this work, it has been assumed that the cell is an 18650 cylindrical cell and therefore 
the ratio Ac/Vc = 253 m'^. Other cell geometries can easily be considered since only the 
convective surface area to cell volume ratio is required in this model. However, for large cells 
the assumption of uniform cell temperature may not be satisfactory as large temperature 
gradients build up within the cell. 

During high C-rate operation, the cell temperature can signihcantly increase and affect 
the cell physical and chemical properties. Therefore, the coupling between the thermal and 
electrochemical model must include the temperature dependency of the model parameters. 
Temperature dependencies of the electrolyte diffusivity and conductivity are taken from [36]. 
The solid phase diffusion coefficient Dg and the reaction kinetics constant k are also highly 
dependent on temperature. A common approach assumes an Arrhenius’ law temperature 
dependency given by [39]: 


^/J 


^""^exp 


R 



(16) 


where ip denotes the parameter considered and ip^^^ is the value of this parameter at . 
Temperature also has an impact on the open-circuit potential of the electrodes. In this 
paper, this was approximated using a hrst-order Taylor series expansion with respect to 
temperature: 

U = U^-f +{T-T^-f) (17) 

where is the open-circuit potential at and {dU / dT) is the entropy change coefficient. 
Empirical expressions reported in [39, 40] for the entropy change as a function of solid-phase 
surface stoichiometry of LiCo02 and MCMB electrodes were used for the simulations. The 
parameters of the thermal model are summarised in Table 2. 


2.3. Chebyshev orthogonal collocation 

The thermal-electrochemical P2D model consists of a set coupled nonlinear partial- 
differential equations in time and space. An analytical solution for such a complex problem 
is not available and numerical methods are employed to spatially discretise the equations in 
the X- and r-directions. The discretised P2D model consists of a system of ODEs and DAEs 
that can be integrated using a standard time-adaptive ODE/DAE solver such as MATLAB’s 
odd5s [44]. The finite difference method has been commonly used to discretise the P2D 
model in space but this requires a significant number of discrete nodes and therefore results 
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Table 2: Thermal model parameters 


Parameter 

Units 

Value 

Ref. 

Cp 

J.kg-\K-^ 

750 

[39] 

P 

kg.m~^ 

1626 

[39] 

rpref 

K 

298 

139] 

h 


30 

- 

Too 

K 

298 

- 


kJ.mol~^ 

35 

[39, 41] 


kJ.mol~^ 

29 

[39, 42] 

Et 

kJ.mol~^ 

20 

[39, 40] 

Et 

kJ.mol~^ 

58 

[39, 43] 


in a large system of equations. In this paper, the electrochemical P2D model is discretised 
using a class of spectral methods called Chebyshev orthogonal collocation that results in a 
much smaller system of equations compared to hnite difference for a similar accuracy [45]. 

Spectral methods consist of expanding the solution m of a differential equation in terms of 
chosen orthogonal basis functions and determining the coefficients of this expansion to satisfy 
the differential equation. For problems with periodic boundary conditions, cosine functions 
are a natural choice of basis functions. However, for non-periodic boundary conditions, 
discontinuities introduced at the boundaries result in Gibbs phenomena that drastically 
impede spectral accuracy. This can be circumvented by adding linear and/or quadratic 
terms to the Fourier series expansion to enforce the boundary conditions as in [23]. However, 
Chebyshev polynomials are a more natural choice for the solution of differential equations 
with non-periodic boundary conditions such as the P2D model [46, 45]. The solution u{x,t) 
of the PDF is therefore approximated by the truncated Chebyshev expansion; 

N 

UN{x,t) = '^Uk{t)Tk{x), a;e[-l,l] (18) 

fc =0 

where Uk{t) are the iV-|-1 Chebyshev coefficients of the expansion that need to be determined 
and Tk denotes the Chebyshev polynomial of the first-kind of degree k. In the present 
work, the coefficients are determined by the so-called orthogonal collocation method, also 
sometimes referred to as the pseudo-spectral method. The coefficients itkit) are calculated 
by forcing the truncated Chebyshev series (18) to satisfy the differential equation exactly at 
the discretising nodes x* given by: 


= cos 1^—J i = 0, l,...,iV (19) 

By choosing the coefficients itkit) so that UNixi,t) = uixi, t), the Chebyshev series expansion 
(18) becomes an interpolating polynomial of degree N to the solution u of the differential 
equation at the collocation nodes Xi. It can be shown [47] that the interpolating polynomial 
can be expressed in terms of the value of the solution at the collocation nodes Ujit) = 



u{xj,t) = Ujq{xj,t) by: 


N 

Un{x, t) = ^ Uj(t)(l)j{x), X e [-1,1] 

j=0 


where the functions 0^ are given by: 




{-iy+'(i-xAT’Ax) 

CjN^ix - Xj) 


X e [-1,1] 


( 20 ) 


( 21 ) 


with Cj = 2 for j = 0 and j = N and Cj = 1 otherwise. When implementing the orthogonal 
collocation method for the solution of PDEs, the coefficients Uk of the series expansion 
are rarely computed explicitly but the differentiation of u is usually performed using a 
differentiation matrix. The pth derivative of the solution u with respect to x evaluated at 
the collocation points can be expressed by: 


u 


ip) I 
N 


N 

i=o 




( 22 ) 


where the coefficients can be determined by evaluating the pth derivative of the in- 

terpolant (20) at the collocation nodes (19). The coefficients are the elements of the 
so-called differentiation matrix D^, which is the discrete approximation to the pth derivative 
operator d^/dx'^. The derivative of u evaluated at the collocation nodes can be expressed 
in terms of the value of u at the collocation nodes u with: 


u(p) = (23) 

The MATLAB function chebdif.m discussed in [47] was used to compute the Chebyshev 
orthogonal collocation differentiation matrices. 

The P2D model must satisfy the boundary conditions associated with the PDEs as dis¬ 
cussed in Section 2.1. In the present work, these boundary conditions are accounted for by 
reducing the size of the differentiation matrix, since each boundary condition gives an addi¬ 
tional constraint that can be used to express the value of the solution at a chosen collocation 
node in terms of the solution values at all the other collocation nodes. This reduces the size 
of the differentiation matrix by one row and one column for each boundary condition consid¬ 
ered. This leads to reduced differentiation matrices that automatically satisfy the boundary 
conditions. 


2.4- Domain decomposition and scaling 

The main advantage of spectral methods over EDM is their fast rate of convergence [45], 
which means that the same accuracy can be obtained with fewer discretisation nodes (re¬ 
duced by a factor of 10 to 100). However, this rapid convergence behaviour, referred to as 
spectral accuracy, is achieved provided that the solution is sufficiently smooth. Discontinu¬ 
ities reduce spectral accuracy and appear in the electrochemical P2D model at each of the 
electrode/separator interfaces. In order to avoid these discontinuities, the cell domain was 
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decomposed into three sub-domains denoted and f2c for the anode, separator and 

cathode sub-domains respectively, where the model equations are solved on distinct sets of 
Chebyshev collocation nodes Xa, Xg and Xc (Fig. 1). Additional interface boundary conditions 
are required to ensure the continuity of the dependent variables and the conservation of flux 
at the interfaces between sub-domains and are summarised in B. The anode and cathode 
solid-phase particles domains flp^a and flp^c respectively are discretised using the same set 
of Chebyshev collocation nodes Xp. The number of Chebyshev collocation nodes used in 
the anode, separator, cathode and particles sub-domains are denoted Na, Ng, Nc, and Np 
respectively. Each of these sub-domains are rescaled to the interval [—1,1] G M, since the 
Chebyshev collocation nodes are dehned on this interval. 

2.5. State-space representation of the discretised model 

The P2D model discretrised by orthogonal collocation consists of a set of non-linear 
differential-algebraic equations (DAEs) with respect to time. Using a state-space represen¬ 
tation, the model can be conveniently written as a semi-explicit DAE system, consisting of 
a set of differential (24) and algebraic equations (25): 

X = f (x, z, u) (24) 

0 = g(x,z,M) (25) 

where the functions f and g are non-linear mapping functions derived from the discretised 
model equations. The state vector x G M”"" associated with the differential equations con¬ 
tains the value at the collocation points of the solid-phase concentration Cg = rcg and the 
electrolyte concentration Cg, as well as the bulk temperature T: 

x=[c^,Ce,T]^ (26) 

The state vector z G associated with the algebraic equations contains the value of the 

volumetric reaction rate at the collocation points and the solid-phase electric potential 
at the cathode and anode current collector 0°^ and 0°^ respectively. 

z = (27) 

The measurement vector y = [U T]^ containing the value of the voltage V and the tem¬ 
perature T is computed from the differential and algebraic state vectors according to the 
measurement equation (28). The input m is a scalar equal to the applied current I. 

H^u (28) 

The derivation of matrices H^, and is trivial since the temperature is a differential 
state of the model and the voltage computation is straightforward from the algebraic state 
vector and the input vector. 

Equations (24), (25) and (28) constitute a state-space representation of the thermal- 
electrochemical P2D model. This representation is particularly convenient from a control 
engineering perspective and can be implemented in the ODE/DAEs MATLAB solver odel5s. 
In addition, this representation is useful for the design and implementation of a state esti¬ 
mator as discussed in Section 4. 


y = [H. H.] 
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Figure 2: Cell voltage under constant-current discharge at several C-rates. Solid lines: COMSOL, markers: 
Chebyshev orthogonal collocation in MATLAB. 

3. Thermal-electrochemical model simulation results and discussion 

In this section, we compare the model prediction obtained from the solution of the P2D 
model solved using the Chebyshev orthogonal collocation method discussed in Section 2.3 to 
the solution obtained using the commercial finite-element software COMSOL Multiphysics. 
The implementation of the thermal-electrochemical P2D model in COMSOL was performed 
using the equation-based modelling toolbox similar to [48], see C. The P2D model solved 
by finite-elements in COMSOL is subsequently referred to as the ‘high-fidelity’model for 
simplicity. 

Fig. 2 compares the cell terminal voltage predicted by both approaches during constant 
current discharge at several C-rates ranging from 1C to IOC. The chosen number of colloca¬ 
tion nodes in the anode, separator and cathode are Na = Q, Ng = 3 and Nc = 6 respectively 
and the number of collocation nodes in each particles of both electrodes is Np = 15. The 
voltage predicted by our approach is in very good agreement with the high-hdelity model up 
to high C-rates (IOC) with a root-mean square and a maximum error of 10 mV and 50 mV 
respectively. The solution of the model using Chebyshev orthogonal collocation is typically 
30 times faster than the solution using COMSOL. The computation of a single discharge 
curve on a desktop computer using a 3.40 GHz processor with 8 GB RAM is performed in 
about 5 min with COMSOL, compared to 1 s to 10 s with our implementation. The order 
of magnitude of these computation times are consistent with results reported in [23] with a 
Maple solver. 

The number of collocation nodes required to discretise the cell domain depends on the 
C-rate, since higher C-rates result in larger gradients of dependent variables across the cell. 
In particular, the accuracy of results highly depends on the number of collocation nodes Np 
in the solid-phase particles. This is due to the very sharp gradients of lithium concentration 
at the surface of these particles for medium to high C-rates. The root-mean square and 
maximum absolute errors between the orthogonal collocation and the high-fidelity model for 
the 1C, 2C and 5C full constant-current discharge with increasing number of nodes in the 
solid-phase particles are shown in Fig. 3a and Fig. 3b respectively. These graphs confirm that 
a larger number of collocation nodes results in smaller error on voltage prediction, and the 
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Figure 3: RMS error (a) and maximum absolute error (b) on voltage, and maximum relative error on bulk 
SOC in both electrodes (c) predicted by the P2D model solved using orthogonal collocation in MATLAB 
compared to the high-fidelity COMSOL model under constant-current discharge at several C-rates with 
respect to the number of collocation nodes in each solid-phase particles Np. 


higher the C-rate the more nodes are required. Although, it is suggested by the maximum 
error graph (Fig. 3b) that more collocation nodes are required for the 1C discharge cycle 
compared to higher C-rate, this is not representative of the whole discharge curve. The 
voltage maximum error arises from the very low SOC portion of the discharge curve, when 
the cell voltage rapidly drops due to the low concentration in the anode material. The 
maximum error tends to be smaller at higher C-rates compared to 1C because such a low 
anode concentration cannot be reached at higher C-rates. 

An important state of the model for a BMS is the cell SOC. In this paper, the bulk SOC 
of the electrode i is dehned according to: 


soa (t) 


^nt) - 

Ql00% _ q0% 


(29) 


where and 6*°^° denote the electrode stoichiometry at 100 % and 0 % respectively. The 
average electrode stoichiometry 9^^^ is calculated by integrating the solid-phase concentration 
in each particle and across the cell according to (30). As shown in Fig. 3c, the maximum 
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relative error on the bulk SOC in both the anode and the cathode rapidly falls below 1 % 
error with less than 10 nodes in the solid-phase particles and below 0.1 % with only 15 nodes 
up to IOC. 


or r 





Cs,i {x, r, t) 

^max 

^s,i 


drdx 


(30) 


These results conhrm the accuracy of our approach. Unlike simpler models that have 
previously been used for battery state estimation, such as the single-particle model, the P2D 
model is able to predict local variations of internal states across the cell. Such variations 
become particularly acute for high C-rate operation (Fig. 4), such as discharge where a 
large amount of lithium is released into the electrolyte at the anode and absorbed at the 
cathode into the active material. Due to the relatively slow diffusion of lithium ions from 
the anode to the cathode, a large concentration gradient builds up in the electrolyte across 
the cell and reduces the cell performance. It can be seen from Fig. 4a that during a IOC 
discharge, the electrolyte is almost depleted of lithium in less than 50 s. The single-particle 
model also assumes that the insertion reaction rate is uniform within each electrode. Fig. 4b 
shows that this assumption is not valid at high C-rates. At short timescales, the reaction 
rate at the electrode-separator interface can be an order of magnitude higher compared 
to the reaction rate at the current collector. The use of the P2D model for battery state 
estimation could provide valuable information on local internal states to the BMS and enable 
the implementation of better health-conscious battery management algorithms [30]. 

In embedded application for automotive BMSs, the state estimation would have to be 
performed with the dynamic current input experienced by the battery pack. In this work, 
we used the Combined ARTEMIS Driving Cycle (CADC) [49] to generate a dynamic current 
excitation prohle approximately representative of an electric car (or PHEV with all-electric 
mode) drive cycle. We assumed that the cell input current was proportional to the vehicle’s 
acceleration and that 25 % of the braking acceleration was recovered to charge the battery. 
We chose the scaling factor between car acceleration and input current in order to obtain 
a relatively aggressive load prohle with peak current reaching 15C. The model prediction 
of voltage and temperature for the CADC input current are shown in Fig. 5. The model 
predicts that the full discharge of the battery occurs in 1700 s and that the temperature would 
rise up to 72 °C under these relatively aggressive conditions. The temperature elevation 
predicted by the model is relatively high compared to what would be experienced by cells in 
an automotive battery pack because of the high peaks of current and the simplistic air cooling 
system considered in this study. However, this simulation demonstrates that the model can 
be solved under highly dynamic and high C-rate operation. The thermal boundary condition 
could easily be changed to replicate a liquid cooling system if required. As illustrated by 
Fig. 6, the main contribution to the global heat generation rate arises from the contact 
resistance heat generation followed by the ohmic heat generation qohm, the reaction heat 
generation qrxn and the reversible heat generation qrev Reversible heat is often neglected in 
the literature for high current operation due to its relatively low magnitude in comparison 
to other heat sources. However, it can be observed that the reversible heat cumulated over 
the driving cycle is not negligible. 

The numerical solution of the model for a dynamic input with high amplitude current 
peaks is more intensive than a constant current discharge, and the computation time required 
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Figure 4: Electrolyte concentration (a) and local volumetric reaction rate (b) profiles computed by the 
thermal-P2D model solved with orthogonal collocation at several time steps estimated under a IOC constant- 
current discharge. 
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Figure 5: Input current (C-rate), voltage and temperature predicted by the thermal-P2D model under a 
Combined ARTEMIS Driving Cycle. 

by the solver is higher compared to constant-current simulations. The solution of 1700 s 
of simulation under the CADC considered required 285 s of computation on the desktop 
computer previously mentioned. However, this is still a relatively low computation time 
since only 168 ms were required on average to solve 1 s of simulation. This is a promising 
result for future work on the real-time solution of the thermal-electrochemical model for 
battery state estimation. 

4. State estimation using a modified EKF 

We now discuss the implementation of a Kalman hlter for the estimation of the model 
states (26) and (27) from noisy measurements of V, I and T. The Kalman hlter is a 
recursive algorithm that infers the battery internal states by correcting the model states in 
order to minimise the error between the predicted voltage and temperature and the actual 
measurements of voltage and temperature for a given current input. 

We hrst give a very brief overview of Kalman hltering for linear models, followed by 
a detailed explanation of the modihcations required for its application to the thermal- 
electrochemical P2D model. These modihcations are motivated by the fact that the P2D 
model (i) is non-linear and (ii) contains algebraic constraints that cannot be handled by 
standard algorithms. In this paper, we applied the extended Kalman hlter algorithm for 
DAEs discussed in [50] to the thermal-electrochemical P2D model. The derivation of the 
Kalman hlter equations is not provided in this paper and we refer the reader instead to the 
literature on Kalman hltering such as [51]. 


15 



















Figure 6: Cumulative heat generated and convective heat removed per cell unit volume under the CADC. 


4-.1. The Kalman filter 

The Kalman filter is a computationally efficient recursive algorithm for state estimation 
of dynamic systems described by a stochastic linear state-space models: 


x(t) = ylx(t) -|- Bu(t) + w(t) (31) 

Yk = Cxfc -h Vfc (32) 


where (31) is the continuous state equation governing the system dynamics, and (32) is 
the discrete measurement equation relating the states of the system x(t) to the available 
measurements y^. The model input u(t) is usually assumed deterministic, while the states 
and measurements are affected by additive uncorrelated zero-mean Gaussian process noise 
w(t) and measurement noises with covariance matrices Q and R respectively, to account 
for random environment disturbances and sensor noise. 

At every time step tk = kTg, with sampling period Tg, an estimate of the state vector x^. 
and the associated error covariance matrix Pk are computed in two steps: the time update 
and the measurement update. In the time update, a priori estimates of the state x^^^ 
and error covariance at tk+i are calculated using the model and the known input u{t) 
according to (33) and (34), where <h = exp {A.Tg) is the state-transition matrix. 

x^+i = $Xfc -h f ^Bu (r) dr (33) 

Jtk 

Pk+i ~ ^Pk + Pk^^ + Q (34) 

In the measurement update, a posteriori estimates of the state x^+i and error covariance 
Pk+i are computed based on the error between estimated measurements and actual 
noisy measurements y^+i according to (35), (36) and (37). 

/W = h+iC-’-(c-Pj;iC-^ + ii)^‘ (35) 

Xfc +1 = + Kk +1 (yfc+i - yfc+i) (36) 

A+i = (/ - Kk+iC) Pfi^, (37) 
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The Kalman filter is the optimal state estimator in the least-squares sense for minimising the 
state estimation error for linear systems. However, the battery model discussed in Section 2 
is non-linear and has algebraic constraints. We therefore discuss a modihed version of the 
Kalman hlter based on the EKF algorithm for the state estimation of the battery model. 

4 . 2 . Battery stochastic state-space model 

The EKF algorithm relies on a non-linear stochastic state-space model. Such a state- 
space representation of the P2D model can be derived from the state-space representation 
given by (24), (25) and (28) by adding process noise and measurement noise to the dynamics 
and measurement equations respectively according to: 

x(t) = f (x(t), z(f), u(t)) + w(t) 

0 = g(x(t),z(t),'u(t)) 

yk = [H^ + HuUk + Vfc 

_Zfc_ 

Similarly to the Kalman hlter discussed in Section 4.1, the process noise w and measure¬ 
ment noise are zero-mean Gaussian additive noises uncorrelated in time with covariance 
matrices Q and R respectively. We assumed no noise on the input current u(t). 

4-3. DAE State-space model linearisation 

The difference between the Kalman hlter and the EKF consists of additional linearisa¬ 
tion steps required between the time update and the measurement update compared to the 
Kalman hlter algorithm for linear models. The linearisation of the diherential equation (38) 
and the algebraic equation (39) are performed about the current state estimate [x, z]^ at 
every time step [50]. This allows the system of non-linear DAEs to be transformed into a 
system of locally linear ODEs that can be used in both the time update and measurement 
update steps discussed in Section 4.1. 

To perform the model linearisation, we hrst dehne the following variables: 


x = x —X (41) 
i = X - i (42) 
z = z — z (43) 
u = u — u (44) 
y = y - y (45) 


Assuming that the functions f and g are sufficiently diherentiable, hrst-order Taylor series 
expansions of these functions about the current state estimate are: 

f (x, z) = f (x, z) -h frX f^Z -h f„h (46) 

g (x, z) = g (x, z)-h g^x-h g^z-h g„h (47) 

where f* and gj denote the partial derivative of f and g respectively with respect to the 
variable i = {x, z,m} evaluated at the current state estimate [x, z]^. The matrices fa,, f2, ga. 


(38) 

(39) 

(40) 
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and gz are therefore Jacobian matrices of the functions f and g. A linear approximation of 
the DAE system is obtained by substituting (46) and (47) into (38) and (39): 


frX -h f^z -h fuU 

(48) 


(49) 

gz;X -F g^Z -h guU 

(50) 


The random variables w and v^. are additive and can therefore be ignored in the linearisation 
process. 

By assuming that the Jacobian matrix is non-singular, which is equivalent to assuming 
that the semi-explicit DAEs system is of index 1, the linearized algebraic constraint (47) can 
be rearranged to obtain an expression of the algebraic state vector in terms of the differential 
state vector: 

z = [gx* + guu] (51) 

Substituting (51) into (49) gives the following linearised state equation that includes the 
algebraic constraint, 

i = (52) 

where: 

= 4 - f.gjig, (53) 

- f.gjig, (54) 


In the linearisation process, the DAE system is therefore transformed into an ODE system 
that can be used in a standard Kalman hlter algorithm. The state-transition matrix <h of the 
linearised model is given by <h = exp and is used in the time update of the Kalman 

hlter discussed in Section 4.1. In a similar way, substituting (51) into the measurement 
equation (40) gives. 


where: 


y = C'*”* + 

(55) 

C-'- = 11,- ffzgf^g. 

(56) 

= H^- Hgf^g^ 

(57) 


The measurement matrix therefore includes the linearised algebraic constraint and can 
be used in the measurement update of the standard Kalman hlter algorithm to compute the 
Kalman gain and update the error covariance estimate. 


4-4- Summary of modified EKF for non-linear DAEs 

This section provides a step-by-step description of the modihed EKF for systems of 
DAEs. In [50], the modihed EKF algorithm is based on the square-root form of the EKF 
for numerical stability. This guarantees that the error covariance matrix remains positive 
semi-dehnite by using the square-root of the error covariance matrix instead of the error 
covariance matrix itself. Although the square-root implementation is more robust in some 
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cases, the two algorithms are mathematically equivalent and we found no difference in the 
results between the standard and square-root form for our problem. Only the standard 
version of the modihed EKF is discussed in this paper for simplicity. We therefore refer the 
reader to [50] regarding the implementation of the square-root EKF. 

The EKF algorithm is initialised by assuming an initial differential state estimate xq 
and error covariance matrix Pq at time to- A consistent initial algebraic state vector zq 
is computed using the MATLAB solver fsolve for systems of non-linear equations. The 
computation of an initial error covariance matrix Pq that accounts for the spatial correlation 
of the model states is crucial for the performance and convergence of the EKF for DAEs. 
The structure of the matrix Pq must satisfy the spatial correlation of the error, otherwise the 
conservation of lithium cannot be guaranteed in the measurement update step. For instance, 
the use of a diagonal initial error covariance matrix (i.e. no spatial correlation) results in 
the DAE solver failure after a few time steps due to violation of conservation laws. In this 
work, the structure of the error covariance matrix was obtained numerically using the model 
states computed by integrating the model under various inputs and initial conditions. 

Once the modihed EKF for DAEs is initialised, the following algorithm steps are per¬ 
formed recursively: 

1. State time update: The current state estimate [x^, z^j is projected forward in time to the 
next time step by integrating the non-linear DAEs system using the MATLAB solver for 
DAEs odelSs from tk to tk+i- The predicted state vector at time tk+i is the a priori state 
estimate \^k+i^^k+-i\^ time ffc+i- 

2. Model linearisation: The DAE model is linearised about the current state estimate [x^, z^j 

to compute the state transition matrix $ = exp of the linearised model. 

3. Error covariance time update: The error covariance Pk is propagated in time using (34) 
to obtain the a priori error covariance estimate Pk+i at time tk+i- 

4. Model linearisation: The DAE model is linearised about the a priori state estimate 
\^k+ii^k+-\^ order to compute the measurement matrix for the measurement 
update. 

5. Measurement update: The matrix previously computed is used to calculate the 

Kalman gain according to (35). The a priori differential state estimate x^+^ and er¬ 
ror covariance estimate are updated to account for the measurement y^+i according 
to (36) and (37) respectively. The measurement estimate y^+i in (36) is computed from 
the prior estimate and the input Uk+i according to the measurement equa¬ 

tion (28). The a posteriori differential state estimate x^+i and error covariance estimate 
Pfc+i are therefore obtained. 

6. Consistent algebraic states The consistent a posteriori algebraic state estimate z^+i is 
obtained from the posterior differential state estimate x^+i and the input Uk+i using the 
MATLAB fsolve function. 

This algorithm is repeated recursively at every time step. 

5. State estimation results and discussion 

In this section, we present simulation results showing the performance of the modi¬ 
hed EKF algorithm discussed in Section 4.4 for the state estimation of the full thermal- 
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Figure 7: Evolution of the voltage (a) and temperature (b) computed by the EKF compared to the actual 
and noisy measurements during the first 200 s of a 4C constant-current discharge. Evolution of the absolute 
error on voltage (c) and anode bulk SOC (d) estimated by the EKF compared to actual values generated by 
the reference simulation during the full 4C constant-current discharge. 

electrochemical model discussed in Section 2. The state estimation of a battery cell requires 
experimental data as inputs to the EKF, namely the applied current and the measured volt¬ 
age and temperature response of the cell. Although the state estimator performance should 
ultimately be tested against real experimental data, due to the difficulty in verifying in situ 
the internal states in a real battery, we used the thermal-electrochemical model itself to 
emulate experimental results. Employing such numerical experiments is worthwhile since 
the state estimate error can be easily computed. 

As a reference case, numerical experiments were computed by integrating the thermal- 
electrochemical model from 100 % SOC until the 2 V minimum cut-off voltage under 
constant-current discharge and the CADC discussed in Section 3. The EKF was then started 
from several initial conditions to check its convergence behaviour at different SOC ranging 
from 100 % to 50 %. For both the CADC and const ant-current tests, the EKF initial guess 
on states assumed a cell at equilibrium (i.e. no concentration gradients) with an error on 
both the anode and cathode SOC of 30 %. The initial error on the temperature was set to 
10°C. The variances for the generation of the additive Gaussian white measurement noise 
were set to cxy = 1 x lO”"^ and = 0.25 for the voltage and temperature respec¬ 
tively. These variances correspond to a standard deviation ay = 10 mV on the voltage and 
aT = 0.5 K on the temperature. The measurement noise covariance matrix R of the EKF 
was dehned using these values. No process noise was added to the state variables of the 
model and therefore the EKF process noise covariance matrix Q was set to zero. 

Fig. 7 shows the voltage and temperature calculated by the EKF and the corresponding 
measurements for the hrst 200 s of a 4C constant-current discharge with a 5 s time step. The 
EKF voltage and temperature rapidly converge to the actual reference values in only a few 
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Figure 8: Evolution of the anode bulk SOC estimated by the EKF compared to the actual value for CADC 
charge/discharge cycle (a). Evolution of the voltage (b), temperature (c) and electrolyte lithium concentra¬ 
tion at the anode and cathode current collectors (d) computed by the EKF compared to the actual and noisy 
values generated by the reference simulation during the first 200 s of the CADC charge/discharge cycle. 


time steps. The voltage absolute error for the full 4C discharge is shown on Fig. 7c. The EKF 
was started at three different initial times, 15 s, 250 s and 500 s, to check the convergence 
behaviour at different SOCs. The grey line represents the 95 % conhdence interval on the 
voltage noisy measurement {2av) and the grey dots are the absolute measurement error on 
voltage. The EKF shows similar convergence behaviour for all initial time to studied and 
the voltage estimate falls below the 95 % conhdence interval within the hrst few time steps. 
Similar results were observed for temperature measurements. 

The EKF algorithm is designed to accurately ht the measurements but this does not 
guarantee the convergence of the state estimates. Fig. 7d shows the absolute error between 
the EKF estimated anode bulk SOC and the actual SOC for the full 4C constant-current 
discharge. Again, the EKF was started at different SOCs (initial times to) and shows a 
satisfactory convergence behaviour for all cases. From the 30 % initial absolute error, the 
anode bulk SOC estimate error falls below 1 % after less than 200 s of simulation (300 s 
when started at to = 250 s). Similar results were observed for the cathode bulk SOC due to 
the conservation of lithium in the cell. 

The EKF algorithm was then applied to the battery charge/discharge cycle under the 
CADC. This drive cycle is highly dynamic with large current peaks. However, the EKF was 
solved using a 5 s time-step to reduce the computation time. The EKF under CADC was 
solved in 0.6 s of computation per second of simulation on average on the desktop computer 
previously mentioned. Fig. 8b and Fig. 8c show the rapid convergence of both the EKF 
voltage and temperature compared to the measurements. Similarly to the constant-current 
scenario, the anode bulk SOC also converges relatively quickly to the actual value. The 
absolute error on anode bulk SOC falls below 1 % by 150 s of simulation. 
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Another interesting state of the model is electrolyte concentration, since saturation or 
depletion of lithium in the electrolyte can lead to battery performance limitations at high 
current peaks. Fig. 8d shows the evolution of the estimated and actual electrolyte con¬ 
centration at the anode and cathode current collectors during the hrst few seconds of the 
driving cycle. This graphs shows that the algorithm is able to recover from the wrong initial 
conditions and track the electrolyte concentration during battery operation accurately. 

6. Conclusion 

A physics-based thermal-electrochemical spatially-distributed pseudo-2D model for lithium- 
ion batteries, so-called Newman model in the literature, is solved with Chebyshev orthogonal 
collocation in MATLAB. This results in a highly reduced number of states and computa¬ 
tion cost compared to commonly employed hnite-difference or hnite-elements methods while 
maintaining accuracy. Comparative results against a much higher-order model solved in 
COMSOL Multiphysics conhrm that accuracy is preserved up to high C-rates. The rela¬ 
tively low number of states required with this approach enables our implementation of the 
model to be combined to a state observer for the estimation of battery internal states. 

We used the extended Kalman hlter to estimate the states of the pseudo-2D battery 
model due to its relatively low computational cost compared to other observers for non-linear 
models. The extended Kalman hlter is able to estimate the state error by using a time-varying 
linear approximation of the model differential-algebraic equations about the state estimate at 
every time-step. To our knowledge, this work is the hrst attempt at estimating the internal 
states of the fully spatially-distributed pseudo-2D battery model using an extended Kalman 
hlter. Results indicates that our state estimation algorithm is able to quickly (less than 200 s) 
recover the states of the model even with a 30 % error on the initial SOC. This approach 
could be used within advanced battery management embedded systems for accurate battery 
state estimation and coupled to degradation models for health-conscious battery control. 
Further work is being undertaken on investigating the observability and identihability of the 
model states and parameters, and developing a parameter estimation algorithm. 
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Appendix A Change of variable for the spherical particle diffusion model 

The spherical diffusion model for the solid-phase particles is modihed using the change 
of variable Cg = rCg. Assuming constant diffusivity, the diffusion equation becomes: 


dcg ^ d'^Cg 

dt ~ " dr^ 


(58) 
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(59) 

(60) 


with the following boundary conditions: 

dcs _ Cs{r = Rs) _ -Rs .Li 

dr Rs asRDs^ 

(r = 0) = 0 

Appendix B Interface boundary conditions for the domain decomposition 

Continuity of electrolyte electric potential 


Lc) 

0e,s(^ La) 

0e,c(^ Lc) 

(61) 

(62) 

Continuity of the concentration profile 



Ce,a(^ La) 
Lq) 

La) 

Lq) 

(63) 

(64) 

Continuity of lithium ion flux 



X — Lid 


(65) 

oIHu) ^ 

x=L^ 

- myw 

X — E/q 

(66) 


Appendix C COMSOL implementation of the thermal-electrochemical P2D 
model 

The COMSOL implementation of the thermal-electrochemical P2D model discussed in 
Section 2.1 is similar to [48] and involves using the COMSOL PDE Interfaces and the 
ODE and DAE Interfaces for equation-based modelling. The macro-scale cell model is 
described on a ID geometry divided into three regions (anode, separator and cathode) in 
the ^(-direction. The micro-scale particle model is described on a 2D geometry in which 
the diffusion coefficients and electronic conductivity are set to zero in the ^-direction [48]. 
This reduces the 2D geometry to a ID geometry in the ^-direction distributed along the 
cell thickness (^(-direction) that is equivalent to the radial r-direction of the solid-phase 
spherical particles of the P2D model. The equations solved on the ID and 2D geometries are 
coupled by projecting the local reaction rate and the solid-phase surface concentration 
c™’’/ from one geometry to the other using the linear extrusion COMSOL function. The ID 
geometry was discretised using a uniformly spaced mesh with 22 elements, 8 elements and 
21 elements in the anode, separator and cathode domains respectively. The 2D geometries 
of the electrodes were discretised using triangular elements for the core of the particles 
and quadrilateral elements at the surface of the particles. The anode was discretised with 
890 triangular elements and 528 quadrilateral elements and the cathode with 588 triangular 
elements and 504 quadrilateral elements. This results in a COMSOL finite-element model 
with 7,856 degrees of freedom. 
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Symbols 

Ac Cell surface area, 

As Electrode surface area, 

as Electrode specific interfacial area, m~^ 

b Bruggeman coefficient 

Ce Electrolyte concentration, mol.m~^ 

Cp Cell lumped specific heat, J.kg~^.K~^ 

Cs Solid-phase concentration, mol.m~^ 

^max Active material max concentration, mol.m~^ 
csurf Solid-phase surface concentration, mol.m~^ 

Dc Electrolyte diffusivity, w?.s~^ 

Ds Solid-phase diffusivity, rn?.s~^ 

E'^ Activation energy of parameter -ip, kJ.moE^ 

T Faraday’s constant, C.mol~^ 

h Convective heat transfer coefficient, W.m~^.K~^ 
I Current, A 

iQ Exchange current density, A.m~‘^ 

iapp Applied current density, A.m~'^ 

ie Current density in electrolyte, A.m~‘^ 
is Current density in solid-phase, A.m~^ 
Volnmetric reaction rate, A.m~^ 
k Reaction rate constant, .s~^ 

q Heat generation rate per unit volume, W.m~^ 

R Gas constant, J.mol~^.K~^ 

Rc Contact resistance, 

Rs Radius of solid-phase particles, m 

T Temperature, K 

Standard state reference temperatnre, K 
Too Coolant temperatnre, K 

Li-ions transference nnmber 
U Electrode open-circnit potential, V 

V Cell voltage, V 

Vc Cell volume, rrA 
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Greek Symbols 

tta Anodic charge transfer coefficient 

ttc Cathodic charge transfer coefficient 

S Thickness of cell layers, m 

Ce Electrolyte volume fraction 

6/ Inert hller volume fraction 

eg Solid-phase volume fraction 

T] Overpotential, V 

K Electrolyte ionic conductivity, S.m~^ 

9s Surface solid-phase stoichiometry 

0 avg Average solid-phase stoichiometry 

6 ^ Initial solid-phase stoichiometry 

p Cell bulk density, kg.m~^ 

a Solid-phase conductivity, S.m~^ 

(pe Electrolyte electric potential, V 

(ps Solid-phase surface electric potential, V 
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